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Abstract 



The bivariate distribution with exponential conditionals (BEC) is introduced by Arnold 
and Strauss [Bivariate distributions with exponential conditionals, J. Amer. Statist. Assoc. 
83 (1988) 522-527]. This work presents a simple and fast algorithm for simulating random 
variates from this density. 
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1 Introduction 

Arnold and Strauss (1988) introduced a bivariate distribution with exponential conditionals 
(BEC), whose unnormalized probability density function is specified by 

/(x,y) = e-('3^+^S'+^^^^^), a;>0, y > 0, (1) 

for P > 0, 7 > 0, 6 > 0. Distribution theory, methods of estimation, and related specifications 
of joint distributions by conditionals can be found in Arnold and Strauss (1988) and Arnold and 
Strauss (1991). The BEC distribution in particular has received attention in applications such 
as reliability analysis (Nadarajah and Kotz, 2006). 

This paper is concerned with simulating random variates from the BEC densities. Arnold 
and Strauss (1988) actually suggested a rejection method using a product exponential as the 
proposal density. While this method is convenient and efficient for some parameter configurations 
(specifically, when 5 in ([1]) is small), it can be shown that as 5 — > oo, its acceptance rate 
approaches zero. To design more efficient algorithms, we may explore several general approaches 
(Devroye, 1986), e.g., inversion, rejection, and ratio of uniforms. Direct inversion is difficult in 
this case because the inverse distribution function of either X or y is not readily available, 
we therefore consider a rejection method. The ratio of uniforms method is considered, but the 
resulting algorithm is not presented here because it also deteriorates as 5 ^ oo and is clearly 
inferior to the proposed method. 

Section 2 presents the new rejection method and evaluates its performance. With a careful 
choice of the envelope function, we obtain an acceptance rate of at least 70%, while keeping the 
algorithm simple and easy to implement. 
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2 A New Rejection Algorithm 

A convenient rejection method, suggested by Arnold and Strauss (1988), is to sample {X,Y) 
from the (unnormalized) proposal density, or envelope function 

go(x,y) = e-'5^-^^ x>0, y > 0, 

and then accept {X,Y) with probability f{X,Y)/gQ{X,Y) = (^-^f^i^'^ , This may be imple- 
mented as the following Algorithm A. 
Algorithm A. 

Step 1. Draw random variates ui,U2,U'i ~ uniform(0, 1) independently. Compute X = — log(ui), 
Y = - \og{u2). 

step 2. If na < e'^^^ return y/7); otherwise go to Step 1. 

For a general rejection algorithm, its acceptance rate is the ratio of the area under /(x, y) 
over that under the envelope function, which for algorithm A simplifies to 

J J 9o[x,y)dxdy Jq 

It is easy to show that Ra{^) ^ 1 as (5 — > and RAiS) — > as 5 — >■ 00. In other words, for large 
(5, the expected number of trials until a pair {X, Y) is accepted can be unreasonably high. 

An alternative strategy is to first sample X according to its marginal density, and then 
sample Y given X. Let us assume for notational convenience /3 = 1 (a simple scaling gives the 
corresponding result for general /?). The (unnormalized) marginal density of X is given by 

fx{x)=e--il + 6x)-\ (2) 

and the conditional of Y given X is 

That is, Y\X ~ exponential 1)/ [7(1 + (^^)]. To accomplish the more difficult part of sampling 
X according to ([2]), consider the function 



g{x; c) 



{l + 5x)~^ 0<x<c 
e~^(l + 6c)^^ X > c 



where c > is a constant to be determined. Clearly 

fx{x) < g{x; c), a; > 0, 

hence g{x; c) is a legitimate envelope function for all c > 0. Drawing X according to g{x; c) is 
simple, because g{x; c) is a mixture whose two components are both easy to sample via inversion. 
Specifically 

g{x; c) = digi{x; c) + ^2^2 (a:; c), 

where 

di = 6-^ log{l + 5c), gi{x; c) = 6/[{l + 6x)\og{l + 6c)], < x < c; 
d2 = e'^/il + 6c), g2{x; c) = e"^+'', x > c. 
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Both gi{x; c) and g2ix; c) arc normalized densities. If we draw X according to gi with 
probabihty di/{di+d2), and according to (72 with the remaining probability, then X is distributed 
according to g overall. This yields the following algorithm for sampling from the original bivariate 
density. 

Algorithm B. 

Step 0. Compute di = S'^ log(l + Sc) and d2 = e~'=/(l + ^c^}. 

Step 1. Draw random variates uq, ui, U2 ~ uniform(0, 1) independently. 

Step 2. When uq < di/{di+d2), set X = + - 1)/6; if U2 < e'^ go to Step 3, otherwise 
go to Step 1. When uq > di/{di + ^2), set X = c — log(tii); if ti2 < (1 + Sc)/{1 + SX) go 
to Step 3, otherwise go to Step 1. 

Step 3. Draw - uniform(0, 1). Return {X/(3, -log{u3)/[-f{l + 6X)]). 

Note that di and ^2 may be pre-compTited if many random variates with the same parameter 
6 are desired. The acceptance rate of algorithm B is easily obtained as 



A natural question is how to determine c. If 5 is small, say 6 < 1, choosing c = 0, which amounts 
to using an exponential envelope, results in a reasonable acceptance rate. Algorithm B reduces 
to Algorithm C in this case. 
Algorithm C. 

Step 1. Draw random variates ui, U2 ^ uniform(0, 1) independently. 
Step 2. Set X = — log(ui); if U2 < {1 + 6X)~^ go to Step 3, otherwise go to Step 1. 
Step 3. Draw U3 ~ uniform(0, 1). Return {X/p, -\og{u3)/['y{l + SX)]). 
The acceptance rate of Algorithm C is 



which coincides with that of Algorithm A. Algorithm C is therefore unsuitable for large 5. (Note 
that, given their identical acceptance rates. Algorithm C has a slight advantage over Algorithm 
A because Algorithm C uses two uniform variates whereas Algorithm A uses three for every 
rejected sample.) 

In contrast, the following shows, for each c > 0, a positive lower bound of the acceptance 
rate of Algorithm B over the range of S. 

Proposition 1. If 6 > and c > then 



Rb{S; c) 



/ fx{x)dx 
J g{x; c)dx 



J„°°e--(l + fe) 



^dx 



di + d2 




Rb{S; c) > {e" + c-^) 



(3) 
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Table 1: Acceptance rate Rb{^', c) of Algorithm B for various values of 6 and c. 



c 





0.1 


0.2 


0.5 


1 


5 

1.5 


2 


3 


5 


10 


20 


100 





1.00 


.916 


.852 


.723 


.596 


.517 


.461 


.386 


.299 


.201 


.130 


.041 


0.5 


.904 


.859 


.829 


.776 


.736 


.719 


.710 


.704 


.705 


.719 


.741 


.796 


0.7 


.836 


.803 


.781 


.747 


.725 


.718 


.716 


.718 


.726 


.746 


.770 


.822 


1 


.731 


.711 


.700 


.684 


.680 


.682 


.687 


.696 


.712 


.737 


.764 


.819 



Proof. We have 

POD 

Rb{6\ c) = {di + ds)"^ / e-^(l + 6x)-^dx 

Jo 

> (di + d2)~^ [ e-%1 + Sx)-^dx 
Jo 

= e-'di{di + d2r\ 

But di/d2 = e'^d-^il + 5c) log (1 + 6c) > ce", where we have used a simple inequality: (1 + 
x) log(l + x) >x when a; > 0. Thus 

Rb{6; c)>e-^di{di + d2T^ 
> e-^ce^ce" + 1)"^ 
= (e'= + c-i)-^ □ 

For large b we need Algorithm B with a good choice of c. Though it is desirable to choose c 
such that -R_b((5; c) is optimized, this optimization is difficult analytically. Time consumed to 
locate the exact maximizer of Rb{S; c) may well offset the improved acceptance rate, especially 
if S changes frequently. Fortunately, it is observed that, when 0.5 < c < 1, Rb{6', c) is quite 
insensitive to the value of c over the full range of 5. Table 1 gives the acceptance rate Rb{6', c) 
for various values of 5 and c = 0, 0.5, 0.7, 1. Note that Algorithm B does not apply if 5 = 0. 
The column for (5 = is taken as lim^^o Rb{S', c). 

With c = 0.7, Rb{^] c) has an approximate lower bound of 0.716. On the other hand, for 
5 < 1, the acceptance rate of Algorithm C, Rc{6) = Rb{6; 0), is bounded below by 0.596. 
Algorithm C has the advantage of simplicity. In addition it avoids the numerical problems of 
Algorithm B when S is near zero. We recommend Algorithm C when 6 < 1 and Algorithm B 
with c = 0.7 otherwise. 
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